! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
! and the University Corporation for Atmospheric Research (UCAR).
!
! Unless noted otherwise source code is licensed under the BSD license.
! Additional copyright and license information can be found in the LICENSE file
! distributed with this code, or at http://mpas-dev.github.com/license.html
!
module atm_core

   use mpas_derived_types
   use mpas_pool_routines
   use mpas_dmpar
   use mpas_log, only : mpas_log_write, mpas_log_info
   use mpas_io_units, only : mpas_new_unit, mpas_release_unit

   !
   ! Abstract interface for routine used to communicate halos of fields
   ! in a named group
   !
   abstract interface
      subroutine halo_exchange_routine(domain, halo_group, ierr)

         use mpas_derived_types, only : domain_type

         type (domain_type), intent(inout) :: domain
         character(len=*), intent(in) :: halo_group
         integer, intent(out), optional :: ierr

      end subroutine halo_exchange_routine
   end interface

   type (MPAS_Clock_type), pointer :: clock


   contains


   function atm_core_init(domain, startTimeStamp) result(ierr)

      use mpas_timekeeping
      use mpas_kind_types
      use mpas_stream_manager
      use mpas_atm_dimensions, only : mpas_atm_set_dims
      use mpas_atm_diagnostics_manager, only : mpas_atm_diag_setup
      use mpas_atm_threading, only : mpas_atm_threading_init
      use atm_time_integration, only : mpas_atm_dynamics_init
      use mpas_timer, only : mpas_timer_start, mpas_timer_stop
      use mpas_attlist, only : mpas_modify_att
      use mpas_string_utils, only : mpas_string_replace
      use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group

      implicit none

      type (domain_type), intent(inout) :: domain
      character(len=*), intent(out) :: startTimeStamp
      integer :: ierr

      real (kind=RKIND), pointer :: dt
      type (block_type), pointer :: block

      logical, pointer :: config_do_restart

      type (mpas_pool_type), pointer :: state
      type (mpas_pool_type), pointer :: mesh
      type (mpas_pool_type), pointer :: diag
      type (field2DReal), pointer :: u_field, pv_edge_field, ru_field, rw_field
      type (field0DReal), pointer :: Time_field
      character (len=StrKIND), pointer :: xtime
      character (len=StrKIND), pointer :: initial_time1, initial_time2
      type (MPAS_Time_Type) :: startTime

      real (kind=RKIND), pointer :: nominalMinDc
      real (kind=RKIND), pointer :: config_len_disp
      real (kind=RKIND), pointer :: Time

      integer, pointer :: nVertLevels, maxEdges, maxEdges2, num_scalars
      character (len=ShortStrKIND) :: init_stream_name
      real (kind=R8KIND) :: input_start_time, input_stop_time


      ierr = 0


      !
      ! Setup threading
      !
      call mpas_atm_threading_init(domain % blocklist, ierr)
      if ( ierr /= 0 ) then
         call mpas_log_write('Threading setup failed for core '//trim(domain % core % coreName), messageType=MPAS_LOG_CRIT)
      end if


      !
      ! Set up inner dimensions used by arrays in optimized dynamics routines
      !
      call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state)
      call mpas_pool_get_dimension(state, 'nVertLevels', nVertLevels)
      call mpas_pool_get_dimension(state, 'maxEdges', maxEdges)
      call mpas_pool_get_dimension(state, 'maxEdges2', maxEdges2)
      call mpas_pool_get_dimension(state, 'num_scalars', num_scalars)
      call mpas_atm_set_dims(nVertLevels, maxEdges, maxEdges2, num_scalars)

      !
      ! Set "local" clock to point to the clock contained in the domain type
      !
      clock => domain % clock
      mpas_log_info => domain % logInfo

      !
      ! Build halo exchange groups and set method for exchanging halos in a group
      !
      call atm_build_halo_groups(domain, ierr)
      if (ierr /= 0) then
         call mpas_log_write('Failed to build halo exchange groups.', messageType=MPAS_LOG_ERR)
         return
      end if

      call mpas_pool_get_config(domain % blocklist % configs, 'config_do_restart', config_do_restart)
      call mpas_pool_get_config(domain % blocklist % configs, 'config_dt', dt)

      !
      ! By default, the 'invariant' stream has an input_interval of "none", so
      ! the following stream read has no effect. However, in case the 'invariant'
      ! stream is defined with an input_interval of "initial_only", we read
      ! time-invariant fields first.
      !
      call MPAS_stream_mgr_read(domain % streamManager, streamID='invariant', whence=MPAS_STREAM_NEAREST, ierr=ierr)
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='invariant', direction=MPAS_STREAM_INPUT, ierr=ierr)

      !
      ! If this is a restart run, read the restart stream, else read the input
      ! stream.
      ! Regardless of which stream we read for initial conditions, reset the
      ! input alarms for both input and restart before reading any remaining
      ! input streams.
      !
      call mpas_timer_start('read_ICs')
      if (config_do_restart) then
         init_stream_name = 'restart'
      else
         init_stream_name = 'input'
      end if
      call mpas_log_write('Reading initial state from '''//trim(init_stream_name)//''' stream')

      call mpas_dmpar_get_time(input_start_time)
      call MPAS_stream_mgr_read(domain % streamManager, streamID=trim(init_stream_name), ierr=ierr)
      call mpas_dmpar_get_time(input_stop_time)
      call mpas_timer_stop('read_ICs')

      if (ierr /= MPAS_STREAM_MGR_NOERR) then
         call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
         call mpas_log_write('Error reading initial conditions',                                                 messageType=MPAS_LOG_ERR)
         call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
      end if

      call mpas_log_write(' Timing for read of '''//trim(init_stream_name)//''' stream: $r s', &
                          realArgs=(/real(input_stop_time - input_start_time, kind=RKIND)/))

      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='input', direction=MPAS_STREAM_INPUT, ierr=ierr)
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_INPUT, ierr=ierr)
      call mpas_log_write(' ----- done reading initial state -----')


      !
      ! Determine horizontal length scale used by horizontal diffusion and 3-d divergence damping
      !
      call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh)
      nullify(nominalMinDc)
      call mpas_pool_get_array(mesh, 'nominalMinDc', nominalMinDc)

      nullify(config_len_disp)
      call mpas_pool_get_config(domain % blocklist % configs, 'config_len_disp', config_len_disp)

      call mpas_log_write('')

      !
      ! If config_len_disp was specified as a valid value, use that
      !
      if (config_len_disp > 0.0_RKIND) then
         call mpas_log_write('Setting nominalMinDc to $r based on namelist option config_len_disp', realArgs=[config_len_disp])

         !
         ! But if nominalMinDc was available in the input file and is different, print a warning
         !
         if (nominalMinDc > 0.0_RKIND .and. abs(nominalMinDc - config_len_disp) > 1.0e-6_RKIND * config_len_disp) then
            call mpas_log_write('nominalMinDc was read from input file as a positive value ($r) that differs', &
                                realArgs=[nominalMinDc], messageType=MPAS_LOG_WARN)
            call mpas_log_write('from the specified config_len_disp value ($r)', &
                                realArgs=[config_len_disp], messageType=MPAS_LOG_WARN)
         end if

         nominalMinDc = config_len_disp

      !
      ! Otherwise, try to use nominalMinDc
      !
      else
         if (nominalMinDc > 0.0_RKIND) then
            call mpas_log_write('Setting config_len_disp to $r based on nominalMinDc value in input file', realArgs=[nominalMinDc])
            config_len_disp = nominalMinDc
         else
            call mpas_log_write('Both config_len_disp and nominalMinDc are <= 0.0.', messageType=MPAS_LOG_ERR)
            call mpas_log_write('Please either specify config_len_disp in the &nhyd_model namelist group,', &
                                messageType=MPAS_LOG_ERR)
            call mpas_log_write('or use an input file that provides a valid value for the nominalMinDc variable.', &
                                messageType=MPAS_LOG_ERR)
            ierr = 1
            return
         end if
      end if


      !
      ! Read all other inputs
      ! For now we don't do this here to match results with previous code; to match requires 
      ! that we read in SST and seaice fields after the call to atm_mpas_init_block()
      !
!      call MPAS_stream_mgr_read(domain % streamManager, ierr=ierr)
!      call MPAS_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_INPUT, ierr=ierr) 

      if (.not. config_do_restart) then
         block => domain % blocklist
         do while (associated(block))
            call mpas_pool_get_subpool(block % structs, 'state', state)
            call mpas_pool_initialize_time_levels(state)
            block => block % next
         end do
      end if

      !
      ! Read a new data stream for Incremental Analysis Update (IAU), if config_IAU_option /= 'off'    : HA (June-15-2016)
      ! FIXME: should I check xtime for the IAU fields? Maybe not.
      ! Note: Because the 'iau' stream has the 'iau' package attached to it, the MPAS_stream_mgr_read( )
      !       call here will actually try to read the stream only if IAU is being used in the run.
      !
      call MPAS_stream_mgr_read(domain % streamManager, streamID='iau', whence=MPAS_STREAM_NEAREST, ierr=ierr)
      if (ierr /= MPAS_STREAM_MGR_NOERR) then
         call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
         call mpas_log_write('Error reading IAU files',                                                          messageType=MPAS_LOG_ERR)
         call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
      end if
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='iau', ierr=ierr)

      !
      ! Set startTimeStamp based on the start time of the simulation clock
      !
      startTime = mpas_get_clock_time(clock, MPAS_START_TIME, ierr)
      call mpas_get_time(startTime, dateTimeString=startTimeStamp) 

      call exchange_halo_group(domain, 'initialization:u')


      !
      ! Perform basic compatibility checks among the fields that were read and the run-time options that were selected
      !
      call mpas_atm_run_compatibility(domain % dminfo, domain % blocklist, domain % streamManager, ierr)
      if (ierr /= 0) then
          call mpas_log_write('Please correct issues with the model input fields and/or namelist.', messageType=MPAS_LOG_ERR)
          return
      end if


      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'mesh', mesh)
         call mpas_pool_get_subpool(block % structs, 'state', state)

         call atm_mpas_init_block(domain % dminfo, domain % streamManager, block, mesh, dt)

         call mpas_pool_get_array(state, 'xtime', xtime, 1)
         xtime = startTimeStamp

         ! Initialize initial_time in second time level. We need to do this because initial state
         ! is read into time level 1, and if we write output from the set of state arrays that
         ! represent the original time level 2, the initial_time field will be invalid.
         call mpas_pool_get_array(state, 'initial_time', initial_time1, 1)
         call mpas_pool_get_array(state, 'initial_time', initial_time2, 2)
         initial_time2 = initial_time1

         ! Set the units to be cf compliant 'seconds since <cf-timestamp>'
         call mpas_pool_get_field(state, 'Time', Time_field)
         call mpas_modify_att(Time_field % attLists(1) % attList, 'units', &
              'seconds since ' // mpas_string_replace(initial_time1, '_', ' '))

          block => block % next
      end do

      call exchange_halo_group(domain, 'initialization:pv_edge,ru,rw')

      call mpas_atm_diag_setup(domain % streamManager, domain % blocklist % configs, &
                               domain % blocklist % structs, domain % clock, domain % dminfo)

      !
      ! Prepare the dynamics for integration
      !
      call mpas_atm_dynamics_init(domain)

   end function atm_core_init


   subroutine atm_simulation_clock_init(core_clock, configs, ierr)

      use mpas_timekeeping

      implicit none

      type (MPAS_Clock_type), intent(inout) :: core_clock
      type (mpas_pool_type), intent(inout) :: configs
      integer, intent(out) :: ierr

      type (MPAS_Time_Type) :: startTime, stopTime
      type (MPAS_TimeInterval_type) :: runDuration, timeStep
      integer :: local_err
      real (kind=RKIND), pointer :: config_dt
      character (len=StrKIND), pointer :: config_start_time
      character (len=StrKIND), pointer :: config_restart_timestamp_name
      character (len=StrKIND), pointer :: config_run_duration
      character (len=StrKIND), pointer :: config_stop_time
      character (len=StrKIND) :: startTimeStamp
      integer :: iounit


      ierr = 0

      call mpas_pool_get_config(configs, 'config_dt', config_dt)
      call mpas_pool_get_config(configs, 'config_start_time', config_start_time)
      call mpas_pool_get_config(configs, 'config_restart_timestamp_name', config_restart_timestamp_name)
      call mpas_pool_get_config(configs, 'config_run_duration', config_run_duration)
      call mpas_pool_get_config(configs, 'config_stop_time', config_stop_time)

      if(trim(config_start_time) == 'file') then
         call mpas_new_unit(iounit)
         open(iounit,file=trim(config_restart_timestamp_name),form='formatted',status='old')
         read(iounit,*) startTimeStamp
         close(iounit)
         call mpas_release_unit(iounit)
      else
        startTimeStamp = config_start_time
      end if
      call mpas_set_time(curr_time=startTime, dateTimeString=startTimeStamp, ierr=local_err)
      call mpas_set_timeInterval(timeStep, dt=config_dt, ierr=local_err)

      if (trim(config_run_duration) /= "none") then
         call mpas_set_timeInterval(runDuration, timeString=config_run_duration, ierr=local_err)
         call mpas_create_clock(core_clock, startTime=startTime, timeStep=timeStep, runDuration=runDuration, ierr=local_err)

         if (trim(config_stop_time) /= "none") then
            call mpas_set_time(curr_time=stopTime, dateTimeString=config_stop_time, ierr=local_err)
            if(startTime + runduration /= stopTime) then
               call mpas_log_write('config_run_duration and config_stop_time are inconsistent: using config_run_duration.', messageType=MPAS_LOG_WARN)
            end if
         end if
      else if (trim(config_stop_time) /= "none") then
         call mpas_set_time(curr_time=stopTime, dateTimeString=config_stop_time, ierr=local_err)
         call mpas_create_clock(core_clock, startTime=startTime, timeStep=timeStep, stopTime=stopTime, ierr=local_err)
      else
          call mpas_log_write('Neither config_run_duration nor config_stop_time were specified.', messageType=MPAS_LOG_ERR)
          ierr = 1
      end if

      !TODO: set phyics alarms here...
      !....
      !....

   end subroutine atm_simulation_clock_init


   subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
   
      use atm_time_integration
      use mpas_rbf_interpolation
      use mpas_vector_reconstruction
      use mpas_stream_manager
      use mpas_atm_boundaries, only : mpas_atm_setup_bdy_masks
#ifdef DO_PHYSICS
!     use mpas_atmphys_aquaplanet
      use mpas_atmphys_control
      use mpas_atmphys_init
      use mpas_atmphys_manager
#endif

      implicit none
   
      type (dm_info), intent(in) :: dminfo
      type (MPAS_streamManager_type), intent(inout) :: stream_manager
      type (block_type), intent(inout) :: block
      type (mpas_pool_type), intent(inout) :: mesh     !MGD does this need to be a pointer?
      real (kind=RKIND), intent(in) :: dt

      type (mpas_pool_type), pointer :: state
      type (mpas_pool_type), pointer :: diag
      type (mpas_pool_type), pointer :: tend
      type (mpas_pool_type), pointer :: sfc_input
      type (mpas_pool_type), pointer :: diag_physics
      type (mpas_pool_type), pointer :: diag_physics_noahmp
      type (mpas_pool_type), pointer :: ngw_input
      type (mpas_pool_type), pointer :: atm_input
      type (mpas_pool_type), pointer :: output_noahmp

      integer :: iCell,iEdge,iVertex,k
      
      real (kind=RKIND), dimension(:,:), pointer :: u, uReconstructX, uReconstructY, uReconstructZ, uReconstructZonal, uReconstructMeridional
      real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, meshScalingDel4
      real (kind=RKIND), dimension(:), pointer :: areaCell, invAreaCell
      real (kind=RKIND), dimension(:), pointer :: dvEdge, invDvEdge
      real (kind=RKIND), dimension(:), pointer :: dcEdge, invDcEdge
      real (kind=RKIND), dimension(:), pointer :: areaTriangle, invAreaTriangle
      integer, pointer :: nCells_ptr, nEdges_ptr, nVertices_ptr, nVertLevels_ptr, nEdgesSolve
      integer :: nCells, nEdges, nVertices, nVertLevels
      integer :: thread
      character(len=StrKIND), pointer :: mminlu

      integer, pointer :: nThreads
      integer, dimension(:), pointer :: cellThreadStart, cellThreadEnd
      integer, dimension(:), pointer :: cellSolveThreadStart, cellSolveThreadEnd
      integer, dimension(:), pointer :: edgeThreadStart, edgeThreadEnd
      integer, dimension(:), pointer :: edgeSolveThreadStart, edgeSolveThreadEnd
      integer, dimension(:), pointer :: vertexThreadStart, vertexThreadEnd
      integer, dimension(:), pointer :: vertexSolveThreadStart, vertexSolveThreadEnd


      logical, pointer :: config_do_restart, config_do_DAcycling

      call atm_compute_signs(mesh)
   
      call mpas_pool_get_subpool(block % structs, 'diag', diag)
      call mpas_pool_get_subpool(block % structs, 'state', state)

      call mpas_pool_get_subpool(block % structs, 'state', state)

      call mpas_pool_get_config(block % configs, 'config_do_restart', config_do_restart)
      call mpas_pool_get_config(block % configs, 'config_do_DAcycling', config_do_DAcycling)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!! Compute inverses to avoid divides later
      
      call mpas_pool_get_array(mesh, 'areaCell', areaCell)
      call mpas_pool_get_array(mesh, 'invAreaCell', invAreaCell)
      
      call mpas_pool_get_array(mesh, 'dvEdge', dvEdge)
      call mpas_pool_get_array(mesh, 'dcEdge', dcEdge)      
      call mpas_pool_get_array(mesh, 'invDvEdge', invDvEdge)
      call mpas_pool_get_array(mesh, 'invDcEdge', invDcEdge)      
      
      call mpas_pool_get_array(mesh, 'areaTriangle', areaTriangle)
      call mpas_pool_get_array(mesh, 'invAreaTriangle', invAreaTriangle)
      
      call mpas_pool_get_dimension(mesh, 'nCells', nCells_ptr)      
      call mpas_pool_get_dimension(mesh, 'nEdges', nEdges_ptr)
      call mpas_pool_get_dimension(mesh, 'nVertices', nVertices_ptr)

      nCells = nCells_ptr
      nEdges = nEdges_ptr
      nVertices = nVertices_ptr
 
      do iCell=1,nCells 
         invAreaCell(iCell) = 1.0_RKIND / areaCell(iCell)
      end do
        
      do iEdge=1,nEdges 
         invDvEdge(iEdge) = 1.0_RKIND / dvEdge(iEdge)
      end do
      
      do iEdge=1,nEdges 
         invDcEdge(iEdge) = 1.0_RKIND / dcEdge(iEdge)
      end do
      
      do iVertex=1,nVertices
         invAreaTriangle(iVertex) = 1.0_RKIND / areaTriangle(iVertex)
      end do
      
      !!!!! End compute inverses
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      call atm_adv_coef_compression(mesh)

      call atm_couple_coef_3rd_order(mesh, block % configs)
      
      call mpas_pool_get_dimension(state, 'nVertices', nVertices_ptr)
      call mpas_pool_get_dimension(state, 'nVertLevels', nVertLevels_ptr)

      nVertices = nVertices_ptr
      nVertLevels = nVertLevels_ptr

      allocate(ke_vertex(nVertLevels,nVertices+1))  ! ke_vertex is a module variable defined in mpas_atm_time_integration.F
      allocate(ke_edge(nVertLevels,nEdges+1))       ! ke_edge is a module variable defined in mpas_atm_time_integration.F
      !$acc parallel default(present)
      !$acc loop vector
      do k = 1, nVertLevels
         ke_vertex(k,nVertices+1) = 0.0_RKIND
         ke_edge(k,nEdges+1) = 0.0_RKIND
      end do
      !$acc end parallel

      call mpas_pool_get_dimension(block % dimensions, 'nThreads', nThreads)

      call mpas_pool_get_dimension(block % dimensions, 'cellThreadStart', cellThreadStart)
      call mpas_pool_get_dimension(block % dimensions, 'cellThreadEnd', cellThreadEnd)
      call mpas_pool_get_dimension(block % dimensions, 'cellSolveThreadStart', cellSolveThreadStart)
      call mpas_pool_get_dimension(block % dimensions, 'cellSolveThreadEnd', cellSolveThreadEnd)

      call mpas_pool_get_dimension(block % dimensions, 'vertexThreadStart', vertexThreadStart)
      call mpas_pool_get_dimension(block % dimensions, 'vertexThreadEnd', vertexThreadEnd)
      call mpas_pool_get_dimension(block % dimensions, 'vertexSolveThreadStart', vertexSolveThreadStart)
      call mpas_pool_get_dimension(block % dimensions, 'vertexSolveThreadEnd', vertexSolveThreadEnd)

      call mpas_pool_get_dimension(block % dimensions, 'edgeThreadStart', edgeThreadStart)
      call mpas_pool_get_dimension(block % dimensions, 'edgeThreadEnd', edgeThreadEnd)
      call mpas_pool_get_dimension(block % dimensions, 'edgeSolveThreadStart', edgeSolveThreadStart)
      call mpas_pool_get_dimension(block % dimensions, 'edgeSolveThreadEnd', edgeSolveThreadEnd)

!$OMP PARALLEL DO
      do thread=1,nThreads
         if (.not. config_do_restart .or. (config_do_restart .and. config_do_DAcycling)) then
            call atm_init_coupled_diagnostics(state, 1, diag, mesh, block % configs, &
                                              cellThreadStart(thread), cellThreadEnd(thread), &
                                              vertexThreadStart(thread), vertexThreadEnd(thread), &
                                              edgeThreadStart(thread), edgeThreadEnd(thread), &
                                              cellSolveThreadStart(thread), cellSolveThreadEnd(thread), &
                                              vertexSolveThreadStart(thread), vertexSolveThreadEnd(thread), &
                                              edgeSolveThreadStart(thread), edgeSolveThreadEnd(thread))
         end if

         call atm_compute_solve_diagnostics(dt, state, 1, diag, mesh, block % configs, &
                                            cellThreadStart(thread), cellThreadEnd(thread), &
                                            vertexThreadStart(thread), vertexThreadEnd(thread), &
                                            edgeThreadStart(thread), edgeThreadEnd(thread))
      end do
!$OMP END PARALLEL DO

      deallocate(ke_vertex)
      deallocate(ke_edge)

      call mpas_rbf_interp_initialize(mesh)
      call mpas_init_reconstruct(mesh)

      call mpas_pool_get_array(state, 'u', u, 1)
      call mpas_pool_get_array(diag, 'uReconstructX', uReconstructX)
      call mpas_pool_get_array(diag, 'uReconstructY', uReconstructY)
      call mpas_pool_get_array(diag, 'uReconstructZ', uReconstructZ)
      call mpas_pool_get_array(diag, 'uReconstructZonal', uReconstructZonal)
      call mpas_pool_get_array(diag, 'uReconstructMeridional', uReconstructMeridional)
      call mpas_reconstruct(mesh, u,                   &
                            uReconstructX,             &
                            uReconstructY,             &
                            uReconstructZ,             &
                            uReconstructZonal,         &
                            uReconstructMeridional     &
                           )
   
#ifdef DO_PHYSICS
      !proceed with initialization of physics parameterization if moist_physics is set to true:
      call mpas_pool_get_subpool(block % structs, 'sfc_input', sfc_input)

      ! Before calling physics_init, ensure that mminlu contains the name of the land use dataset
      call mpas_pool_get_array(sfc_input, 'mminlu', mminlu)
      if (len_trim(mminlu) == 0) then
            call mpas_log_write('****************************************************************')
            call mpas_log_write('No information on land use dataset is available.')
            call mpas_log_write('Assume that we are using ''USGS''.')
            call mpas_log_write('****************************************************************')
            write(mminlu,'(a)') 'USGS'
      end if


      if (moist_physics) then
         !initialization of some input variables in registry:
         call mpas_pool_get_subpool(block % structs, 'tend', tend)
         call mpas_pool_get_subpool(block % structs, 'diag_physics', diag_physics)
         call mpas_pool_get_subpool(block % structs, 'diag_physics_noahmp', diag_physics_noahmp)
         call mpas_pool_get_subpool(block % structs, 'ngw_input', ngw_input)
         call mpas_pool_get_subpool(block % structs, 'atm_input', atm_input)
         call mpas_pool_get_subpool(block % structs, 'output_noahmp', output_noahmp)
         call physics_tables_init(dminfo, block % configs)
         call physics_registry_init(mesh, block % configs, sfc_input)
         call physics_run_init(block % configs, mesh, state, clock, stream_manager)

         !initialization of all physics:
         call physics_init(dminfo, stream_manager, clock, block % configs, mesh, diag, tend, state, 1,       &
                           diag_physics, diag_physics_noahmp, ngw_input, atm_input, sfc_input, output_noahmp)
      endif
#endif
   
      call atm_compute_mesh_scaling(mesh, block % configs)

      call atm_compute_damping_coefs(mesh, block % configs)

      !
      ! Set up mask fields used in limited-area simulations
      !
      call mpas_atm_setup_bdy_masks(mesh, block % configs)

      call mpas_pool_get_dimension(mesh, 'nEdgesSolve', nEdgesSolve)
      call mpas_pool_get_array(mesh, 'meshScalingDel2', meshScalingDel2)
      call mpas_pool_get_array(mesh, 'meshScalingDel4', meshScalingDel4)

      call mpas_log_write('min/max of meshScalingDel2 = $r $r', &
                          realArgs=(/minval(meshScalingDel2(1:nEdgesSolve)), maxval(meshScalingDel2(1:nEdgesSolve))/))
      call mpas_log_write('min/max of meshScalingDel4 = $r $r', &
                          realArgs=(/minval(meshScalingDel4(1:nEdgesSolve)), maxval(meshScalingDel4(1:nEdgesSolve))/))

   end subroutine atm_mpas_init_block
   
   
   function atm_core_run(domain) result(ierr)
   
      use mpas_timekeeping
      use mpas_kind_types
      use mpas_stream_manager
      use mpas_derived_types, only : MPAS_STREAM_LATEST_BEFORE, MPAS_STREAM_INPUT, MPAS_STREAM_INPUT_OUTPUT
      use mpas_timer, only : mpas_timer_start, mpas_timer_stop
      use mpas_atm_boundaries, only : mpas_atm_update_bdy_tend
      use mpas_atm_diagnostics_manager, only : mpas_atm_diag_update, mpas_atm_diag_compute, mpas_atm_diag_reset
   
      implicit none
   
      type (domain_type), intent(inout) :: domain
      integer :: ierr
   
      real (kind=RKIND) :: dt
      logical, pointer :: config_do_restart
      logical, pointer :: config_apply_lbcs
      type (block_type), pointer :: block_ptr

      type (MPAS_Time_Type) :: currTime
      character(len=StrKIND) :: timeStamp
      character (len=StrKIND), pointer :: config_restart_timestamp_name
      integer :: itimestep
      integer :: iounit

      integer :: stream_dir
      character(len=StrKIND) :: input_stream, read_time

      type (mpas_pool_type), pointer :: state, diag, mesh, diag_physics, tend, tend_physics

      ! For timing information
      real (kind=R8KIND) :: integ_start_time, integ_stop_time
      real (kind=R8KIND) :: diag_start_time, diag_stop_time
      real (kind=R8KIND) :: input_start_time, input_stop_time
      real (kind=R8KIND) :: output_start_time, output_stop_time
      
      ierr = 0

      clock => domain % clock
      mpas_log_info => domain % logInfo
      
      call mpas_get_timeInterval(mpas_get_clock_timestep(clock, ierr), dt=dt)
      call mpas_pool_get_config(domain % blocklist % configs, 'config_do_restart', config_do_restart)
      call mpas_pool_get_config(domain % blocklist % configs, 'config_restart_timestamp_name', config_restart_timestamp_name)

      ! Avoid writing a restart file at the initial time
      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_OUTPUT, ierr=ierr)

      ! Also, for restart runs, avoid writing the initial history or diagnostics fields to avoid overwriting those from the preceding run
      if (config_do_restart) then
         call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='output', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
         call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='diagnostics', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
      end if

      call mpas_dmpar_get_time(diag_start_time)

      if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
         block_ptr => domain % blocklist
         do while (associated(block_ptr))
            call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
            call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
#ifdef DO_PHYSICS
            call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
#endif
            call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
            call atm_compute_output_diagnostics(state, 1, diag, mesh)

            block_ptr => block_ptr % next
         end do
      end if
      call mpas_timer_start('diagnostic_fields')
      call mpas_atm_diag_reset()
      call mpas_atm_diag_update()
      call mpas_atm_diag_compute()
      call mpas_timer_stop('diagnostic_fields')

      call mpas_dmpar_get_time(diag_stop_time)

      call mpas_timer_start('stream_output')
      call mpas_dmpar_get_time(output_start_time)
      call mpas_stream_mgr_write(domain % streamManager, ierr=ierr)
      call mpas_dmpar_get_time(output_stop_time)
      call mpas_timer_stop('stream_output')
      if (ierr /= MPAS_STREAM_MGR_NOERR .and. &
          ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_FILE .and. &
          ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_REC) then
         call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
         call mpas_log_write('Error writing one or more output streams',                                         messageType=MPAS_LOG_ERR)
         call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
      end if
      if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
         call mpas_log_write('Timing for diagnostic computation: $r s', realArgs=(/real(diag_stop_time - diag_start_time, kind=RKIND)/))
         call mpas_log_write('Timing for stream output: $r s', realArgs=(/real(output_stop_time - output_start_time, kind=RKIND)/))
      end if

      call mpas_timer_start('diagnostic_fields')
      call mpas_atm_diag_reset()
      call mpas_timer_stop('diagnostic_fields')

      call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)

      block_ptr => domain % blocklist
      call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
      call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
      call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
#ifdef DO_PHYSICS
      call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
#endif

      call mpas_pool_get_config(domain % blocklist % configs, 'config_apply_lbcs', config_apply_lbcs)

      !
      ! Read initial boundary state
      !
      if (config_apply_lbcs .and. &
          MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)) then
         block_ptr => domain % blocklist
         do while (associated(block_ptr))
            call mpas_atm_update_bdy_tend(clock, domain % streamManager, block_ptr, .true., ierr)
            if (ierr /= 0) then
               currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
               call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp)
               call mpas_log_write('Failed to process LBC data on or before '//trim(timeStamp), messageType=MPAS_LOG_ERR)
               return
            end if

            block_ptr => block_ptr % next
         end do
      end if

      ! During integration, time level 1 stores the model state at the beginning of the
      !   time step, and time level 2 stores the state advanced dt in time by timestep(...)
      itimestep = 1
      do while (.not. mpas_is_clock_stop_time(clock))

         currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
         call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)         

         call mpas_log_write('')
         call mpas_log_write('Begin timestep '//trim(timeStamp))

         !
         ! Read future boundary state and compute boundary tendencies
         !
         if (config_apply_lbcs .and. &
             MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)) then
            block_ptr => domain % blocklist
            do while (associated(block_ptr))
               call mpas_atm_update_bdy_tend(clock, domain % streamManager, block_ptr, .false., ierr)
               if (ierr /= 0) then
                  call mpas_log_write('Failed to process LBC data at next time after '//trim(timeStamp), messageType=MPAS_LOG_ERR)
                  return
               end if

               block_ptr => block_ptr % next
            end do
         end if

         ! Regardless of whether boundary tendencies were updated, above, we do not want to read the 'lbc_in' stream
         ! as a general input stream, below.
         call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='lbc_in', direction=MPAS_STREAM_INPUT, ierr=ierr)


         !
         ! Read external field updates
         !
         call MPAS_stream_mgr_begin_iteration(domain % streamManager, ierr=ierr)
         do while (MPAS_stream_mgr_get_next_stream(domain % streamManager, streamID=input_stream, directionProperty=stream_dir))
            if (stream_dir == MPAS_STREAM_INPUT .or. stream_dir == MPAS_STREAM_INPUT_OUTPUT) then
               if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID=input_stream, &
                                                  direction=MPAS_STREAM_INPUT, ierr=ierr)) then
                  call mpas_timer_start('stream_input')
                  call mpas_dmpar_get_time(input_start_time)
                  call MPAS_stream_mgr_read(domain % streamManager, streamID=input_stream, whence=MPAS_STREAM_LATEST_BEFORE, &
                                            actualWhen=read_time, ierr=ierr)
                  call mpas_dmpar_get_time(input_stop_time)
                  call mpas_timer_stop('stream_input')
                  if (ierr /= MPAS_STREAM_MGR_NOERR) then
                     call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
                     call mpas_log_write('Error reading input stream '//trim(input_stream),                                  messageType=MPAS_LOG_ERR)
                     call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
                  end if

                  call mpas_log_write('----------------------------------------------------------------------')
                  call mpas_log_write('  Read '''//trim(input_stream)//''' input stream valid at '//trim(read_time))
                  call mpas_log_write('   Timing for stream input: $r s', realArgs=(/real(input_stop_time - input_start_time, kind=RKIND)/))
                  call mpas_log_write('----------------------------------------------------------------------')

                  call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID=input_stream, direction=MPAS_STREAM_INPUT, ierr=ierr)
               end if
            end if
         end do

         call mpas_timer_start("time integration")
         call mpas_dmpar_get_time(integ_start_time)
         call atm_do_timestep(domain, dt, itimestep)
         call mpas_dmpar_get_time(integ_stop_time)
         call mpas_timer_stop("time integration")   
         call mpas_log_write(' Timing for integration step: $r s', realArgs=(/real(integ_stop_time - integ_start_time, kind=RKIND)/))

         ! Move time level 2 fields back into time level 1 for next time step
         call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state)
         call mpas_pool_shift_time_levels(state)
         
         ! Advance clock before writing output
         itimestep = itimestep + 1
         call mpas_advance_clock(clock)
         currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)


         call mpas_dmpar_get_time(diag_start_time)

         !
         ! Write any output streams that have alarms ringing, after computing diagnostics fields
         !
         call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
         if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
           block_ptr => domain % blocklist
           do while (associated(block_ptr))
              call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
              call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
#ifdef DO_PHYSICS
              call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
              call mpas_pool_get_subpool(block_ptr % structs, 'tend_physics', tend_physics)
#endif
              call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
              call mpas_pool_get_subpool(block_ptr % structs, 'tend', tend)
              call atm_compute_output_diagnostics(state, 1, diag, mesh)

              block_ptr => block_ptr % next
           end do
         end if

         call mpas_timer_start('diagnostic_fields')
         call mpas_atm_diag_update()
         call mpas_atm_diag_compute()
         call mpas_timer_stop('diagnostic_fields')
         call mpas_dmpar_get_time(diag_stop_time)

         call mpas_timer_start('stream_output')
         call mpas_dmpar_get_time(output_start_time)
         call mpas_stream_mgr_write(domain % streamManager, ierr=ierr)
         call mpas_dmpar_get_time(output_stop_time)
         call mpas_timer_stop('stream_output')
         if (ierr /= MPAS_STREAM_MGR_NOERR .and. &
             ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_FILE .and. &
             ierr /= MPAS_STREAM_MGR_ERR_CLOBBER_REC) then
            call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_ERR)
            call mpas_log_write('Error writing one or more output streams',                                         messageType=MPAS_LOG_ERR)
            call mpas_log_write('********************************************************************************', messageType=MPAS_LOG_CRIT)
         end if
         if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
            call mpas_log_write('Timing for diagnostic computation: $r s', realArgs=(/real(diag_stop_time - diag_start_time, kind=RKIND)/))
            call mpas_log_write('Timing for stream output: $r s', realArgs=(/real(output_stop_time - output_start_time, kind=RKIND)/))
         end if

         ! reset any diagnostics here

         if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='diagnostics', direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
            block_ptr => domain % blocklist
            do while (associated(block_ptr))

#ifdef DO_PHYSICS
               call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
#endif
               call atm_reset_diagnostics(diag_physics)

               block_ptr => block_ptr % next
            end do
         end if


         ! Only after we've successfully written the restart file should we we
         !    write the restart_timestamp file
         if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_OUTPUT, ierr=ierr)) then
            if (domain % dminfo % my_proc_id == 0) then
               call mpas_new_unit(iounit)
               open(iounit,file=trim(config_restart_timestamp_name),form='formatted',status='replace')
               write(iounit,*) trim(timeStamp)
               close(iounit)
               call mpas_release_unit(iounit)
            end if
         end if

         call mpas_timer_start('diagnostic_fields')
         call mpas_atm_diag_reset()
         call mpas_timer_stop('diagnostic_fields')

         call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr)

      end do
   
   end function atm_core_run

   
   subroutine atm_compute_output_diagnostics(state, time_lev, diag, mesh)
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Compute diagnostic fields for a domain to be written to history files
   !
   ! Input: state - contains model prognostic fields
   !        mesh  - contains grid metadata
   !
   ! Output: state - upon returning, diagnostic fields will have be computed
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
      use mpas_constants
   
      implicit none
   
      type (mpas_pool_type), intent(inout) :: state
      integer, intent(in) :: time_lev            ! which time level to use from state
      type (mpas_pool_type), intent(inout) :: diag
      type (mpas_pool_type), intent(in) :: mesh
   
      integer :: iCell, k
      integer, pointer :: nCells, nVertLevels, index_qv
      real (kind=RKIND), dimension(:,:), pointer :: theta, rho, theta_m, rho_zz, zz
      real (kind=RKIND), dimension(:,:), pointer :: pressure_base, pressure_p, pressure
      real (kind=RKIND), dimension(:,:,:), pointer :: scalars

      call mpas_pool_get_dimension(mesh, 'nCells', nCells)
      call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
      call mpas_pool_get_dimension(state, 'index_qv', index_qv)

      call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev)
      call mpas_pool_get_array(state, 'rho_zz', rho_zz, time_lev)
      call mpas_pool_get_array(state, 'scalars', scalars, time_lev)

      call mpas_pool_get_array(diag, 'theta', theta)
      call mpas_pool_get_array(diag, 'rho', rho)
      call mpas_pool_get_array(diag, 'pressure_p', pressure_p)
      call mpas_pool_get_array(diag, 'pressure_base', pressure_base)
      call mpas_pool_get_array(diag, 'pressure', pressure)

      call mpas_pool_get_array(mesh, 'zz', zz)

      do iCell=1,nCells
         do k=1,nVertLevels
            theta(k,iCell) = theta_m(k,iCell) / (1._RKIND + rvord * scalars(index_qv,k,iCell))
            rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell)
            pressure(k,iCell) = pressure_base(k,iCell) + pressure_p(k,iCell)
         end do
      end do

   end subroutine atm_compute_output_diagnostics
   
   
   subroutine atm_reset_diagnostics(diag_physics)
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! reset some diagnostics after output
   !
   ! Input: diag_physics  - contains physics diagnostic fields
   !
   ! Output: whatever diagnostics need resetting after output
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
      implicit none
   
      type (mpas_pool_type), pointer :: diag_physics
   
      real (kind=RKIND), dimension(:), pointer :: refl10cm_1km_max

#ifdef DO_PHYSICS
      call mpas_pool_get_array(diag_physics, 'refl10cm_1km_max', refl10cm_1km_max)
      if(associated(refl10cm_1km_max)) then
         refl10cm_1km_max(:) = 0.
      endif
#endif
   
   end subroutine atm_reset_diagnostics


   subroutine atm_do_timestep(domain, dt, itimestep)
   
      use mpas_timekeeping
      use mpas_kind_types
      use atm_time_integration
#ifdef DO_PHYSICS
      use mpas_atmphys_control
      use mpas_atmphys_driver
      use mpas_atmphys_manager
      use mpas_atmphys_update
#endif
      use mpas_atm_halos, only: exchange_halo_group
   
      implicit none
   
      type (domain_type), intent(inout) :: domain 
      real (kind=RKIND), intent(in) :: dt
      integer, intent(in) :: itimestep
      
      type (MPAS_Time_Type) :: startTime, currTime
      type (MPAS_TimeInterval_Type) :: xtimeTime
      integer :: s, s_n, s_d
      real (kind=RKIND) :: xtime_s
      integer :: ierr

      clock => domain % clock
      mpas_log_info => domain % logInfo

      startTime = mpas_get_clock_time(clock, MPAS_START_TIME, ierr)
      currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
         
      xtimeTime = currTime - startTime
      call mpas_get_timeInterval(interval=xtimeTime, S=s, S_n=s_n, S_d=s_d, ierr=ierr)         
      xtime_s = (s + s_n / s_d)


#ifdef DO_PHYSICS
      !proceed with physics if moist_physics is set to true:
      if(moist_physics) then
         call physics_timetracker(domain,dt,clock,itimestep,xtime_s)
         call physics_driver(domain,itimestep,xtime_s)
      endif
#endif

      call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group)

   end subroutine atm_do_timestep
   
   
   function atm_core_finalize(domain) result(ierr)
   
      use mpas_decomp
      use mpas_timekeeping
      use mpas_atm_diagnostics_manager, only : mpas_atm_diag_cleanup
      use mpas_atm_threading, only : mpas_atm_threading_finalize
      use atm_time_integration, only : mpas_atm_dynamics_finalize
      use mpas_atm_halos, only: atm_destroy_halo_groups
   
#ifdef DO_PHYSICS
      use mpas_atmphys_finalize
#endif

      implicit none
   
      type (domain_type), intent(inout) :: domain
      type (block_type), pointer :: block_ptr 
      integer :: ierr

      ierr = 0

      clock => domain % clock
      mpas_log_info => domain % logInfo

      !
      ! Finalize the dynamics
      !
      call mpas_atm_dynamics_finalize(domain)

      call mpas_atm_diag_cleanup()

      call mpas_destroy_clock(clock, ierr)
      call mpas_decomp_destroy_decomp_list(domain % decompositions)

#ifdef DO_PHYSICS
         block_ptr => domain % blocklist
         do while (associated(block_ptr))
            call atmphys_finalize(block_ptr%configs)

            block_ptr => block_ptr%next
         end do
#endif

      !
      ! Destroy halo exchange groups
      !
      call atm_destroy_halo_groups(domain, ierr)
      if (ierr /= 0) then
         call mpas_log_write('Failed to destroy halo exchange groups.', messageType=MPAS_LOG_ERR)
      end if

      !
      ! Finalize threading
      !
      call mpas_atm_threading_finalize(domain % blocklist)

      call mpas_log_write('')
      call mpas_log_write('********************************************************')
      call mpas_log_write('   Finished running the atmosphere core')
      call mpas_log_write('********************************************************')
   
   end function atm_core_finalize


   subroutine atm_compute_mesh_scaling(mesh, configs)

      implicit none

      type (mpas_pool_type), intent(inout) :: mesh
      type (mpas_pool_type), intent(in) :: configs

      integer :: iCell,iEdge, cell1, cell2
      integer, pointer :: nEdges, nCells
      integer, dimension(:,:), pointer :: cellsOnEdge
      real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
      real (kind=RKIND), dimension(:), pointer :: meshScalingRegionalCell, meshScalingRegionalEdge
      logical, pointer :: config_h_ScaleWithMesh

      call mpas_pool_get_array(mesh, 'meshDensity', meshDensity)
      call mpas_pool_get_array(mesh, 'meshScalingDel2', meshScalingDel2)
      call mpas_pool_get_array(mesh, 'meshScalingDel4', meshScalingDel4)
      call mpas_pool_get_array(mesh, 'meshScalingRegionalCell', meshScalingRegionalCell)
      call mpas_pool_get_array(mesh, 'meshScalingRegionalEdge', meshScalingRegionalEdge)
      call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)

      call mpas_pool_get_dimension(mesh, 'nEdges', nEdges)
      call mpas_pool_get_dimension(mesh, 'nCells', nCells)

      call mpas_pool_get_config(configs, 'config_h_ScaleWithMesh', config_h_ScaleWithMesh)

      !
      ! Compute the scaling factors to be used in the del2 and del4 dissipation
      !
      meshScalingDel2(:) = 1.0
      meshScalingDel4(:) = 1.0
      if (config_h_ScaleWithMesh) then
         do iEdge=1,nEdges
            cell1 = cellsOnEdge(1,iEdge)
            cell2 = cellsOnEdge(2,iEdge)
            meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**0.25
            meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**0.75
         end do
      end if

      !
      ! Compute the scaling factors to be used in relaxation zone of regional configuration
      !
      meshScalingRegionalCell(:) = 1.0
      meshScalingRegionalEdge(:) = 1.0
      if (config_h_ScaleWithMesh) then
         do iEdge=1,nEdges
            cell1 = cellsOnEdge(1,iEdge)
            cell2 = cellsOnEdge(2,iEdge)
            meshScalingRegionalEdge(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**0.25
         end do

         do iCell=1,nCells
            meshScalingRegionalCell(iCell) = 1.0 / (meshDensity(iCell))**0.25
         end do
      end if

   end subroutine atm_compute_mesh_scaling


   subroutine atm_compute_signs(mesh)

      implicit none

      type (mpas_pool_type), intent(inout) :: mesh

      integer :: i, j, iCell, iVtx
      integer, pointer :: nCells, nVertices, nEdges, vertexDegree
      integer, dimension(:), pointer :: nEdgesOnCell
      integer, dimension(:,:), pointer :: edgesOnVertex, verticesOnEdge, cellsOnEdge, edgesOnCell
      integer, dimension(:,:), pointer :: verticesOnCell, cellsOnVertex, kiteForCell
      real (kind=RKIND), dimension(:,:), pointer :: edgesOnVertex_sign, edgesOnCell_sign
      real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3, zb_cell, zb3_cell


      call mpas_pool_get_array(mesh, 'edgesOnVertex', edgesOnVertex)
      call mpas_pool_get_array(mesh, 'verticesOnEdge', verticesOnEdge)
      call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
      call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
      call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
      call mpas_pool_get_array(mesh, 'verticesOnCell', verticesOnCell)
      call mpas_pool_get_array(mesh, 'cellsOnVertex', cellsOnVertex)
      call mpas_pool_get_array(mesh, 'edgesOnVertex_sign', edgesOnVertex_sign)
      call mpas_pool_get_array(mesh, 'edgesOnCell_sign', edgesOnCell_sign)
      call mpas_pool_get_array(mesh, 'kiteForCell', kiteForCell)
      call mpas_pool_get_array(mesh, 'zb', zb)
      call mpas_pool_get_array(mesh, 'zb3', zb3)
      call mpas_pool_get_array(mesh, 'zb_cell', zb_cell)
      call mpas_pool_get_array(mesh, 'zb3_cell', zb3_cell)

      call mpas_pool_get_dimension(mesh, 'nCells', nCells)
      call mpas_pool_get_dimension(mesh, 'nVertices', nVertices)
      call mpas_pool_get_dimension(mesh, 'nEdges', nEdges)
      call mpas_pool_get_dimension(mesh, 'vertexDegree', vertexDegree)


      do iVtx=1,nVertices
         do i=1,vertexDegree
            if (edgesOnVertex(i,iVtx) <= nEdges) then
               if (iVtx == verticesOnEdge(2,edgesOnVertex(i,iVtx))) then
                  edgesOnVertex_sign(i,iVtx) = 1.0
               else
                  edgesOnVertex_sign(i,iVtx) = -1.0
               end if
            else
               edgesOnVertex_sign(i,iVtx) = 0.0
            end if
         end do
      end do

      do iCell=1,nCells
         do i=1,nEdgesOnCell(iCell)
            if (edgesOnCell(i,iCell) <= nEdges) then
               if (iCell == cellsOnEdge(1,edgesOnCell(i,iCell))) then
                  edgesOnCell_sign(i,iCell) = 1.0
                  zb_cell(:,i,iCell) = zb(:,1,edgesOnCell(i,iCell))
                  zb3_cell(:,i,iCell) = zb3(:,1,edgesOnCell(i,iCell))
               else
                  edgesOnCell_sign(i,iCell) = -1.0
                  zb_cell(:,i,iCell) = zb(:,2,edgesOnCell(i,iCell))
                  zb3_cell(:,i,iCell) = zb3(:,2,edgesOnCell(i,iCell))
               end if
            else
               edgesOnCell_sign(i,iCell) = 0.0
            end if
         end do
      end do

      do iCell=1,nCells
         do i=1,nEdgesOnCell(iCell)
            iVtx = verticesOnCell(i,iCell)
            if (iVtx <= nVertices) then
               do j=1,vertexDegree
                  if (iCell == cellsOnVertex(j,iVtx)) then
                     kiteForCell(i,iCell) = j
                     exit
                  end if
               end do
               if (j > vertexDegree) then
                  call mpas_log_write('Unexpected error while identifying kiteForCell', messageType=MPAS_LOG_ERR)
               end if
            else
               kiteForCell(i,iCell) = 1
            end if
         end do
      end do

   end subroutine atm_compute_signs


   subroutine atm_compute_damping_coefs(mesh, configs)

      implicit none

      type (mpas_pool_type), intent(inout) :: mesh
      type (mpas_pool_type), intent(in) :: configs

      integer :: iCell, k
      integer, pointer :: nCells, nVertLevels
      real (kind=RKIND), pointer :: config_xnutr, config_zd
      real (kind=RKIND) :: z, zt, m1, pii
      real (kind=RKIND), dimension(:,:), pointer :: dss, zgrid
      real (kind=RKIND), dimension(:), pointer :: meshDensity
      real (kind=RKIND) :: dx_scale_power

      m1 = -1.0
      pii = acos(m1)

      call mpas_pool_get_dimension(mesh, 'nCells', nCells)
      call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)

      call mpas_pool_get_array(mesh, 'meshDensity', meshDensity)
      call mpas_pool_get_array(mesh, 'dss', dss)
      call mpas_pool_get_array(mesh, 'zgrid', zgrid)

      call mpas_pool_get_config(configs, 'config_zd', config_zd)
      call mpas_pool_get_config(configs, 'config_xnutr', config_xnutr)

      dx_scale_power = 1.0
      dss(:,:) = 0.0
      do iCell=1,nCells
         zt = zgrid(nVertLevels+1,iCell)
         do k=1,nVertLevels
            z = 0.5*(zgrid(k,iCell) + zgrid(k+1,iCell))
            if (z > config_zd) then
               dss(k,iCell) = config_xnutr*sin(0.5*pii*(z-config_zd)/(zt-config_zd))**2.0
               dss(k,iCell) = dss(k,iCell) / meshDensity(iCell)**(0.25*dx_scale_power)
            end if
         end do
      end do

   end subroutine atm_compute_damping_coefs


   subroutine atm_adv_coef_compression( mesh )

      implicit none

      type (mpas_pool_type), intent(inout) :: mesh


      real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
      real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
      integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge
      integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge
      real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge

      integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
      integer, pointer :: nCells, nEdges
      integer :: cell_list(20)
      logical :: addcell


      call mpas_pool_get_array(mesh, 'deriv_two', deriv_two)
      call mpas_pool_get_array(mesh, 'adv_coefs', adv_coefs)
      call mpas_pool_get_array(mesh, 'adv_coefs_3rd', adv_coefs_3rd)
      call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
      call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
      call mpas_pool_get_array(mesh, 'advCellsForEdge', advCellsForEdge)
      call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
      call mpas_pool_get_array(mesh, 'nAdvCellsForEdge', nAdvCellsForEdge)
      call mpas_pool_get_array(mesh, 'dcEdge', dcEdge)
      call mpas_pool_get_array(mesh, 'dvEdge', dvEdge)

      call mpas_pool_get_dimension(mesh, 'nCells', nCells)
      call mpas_pool_get_dimension(mesh, 'nEdges', nEdges)

      do iEdge = 1, nEdges
         nAdvCellsForEdge(iEdge) = 0
         cell1 = cellsOnEdge(1,iEdge)
         cell2 = cellsOnEdge(2,iEdge)
         !
         ! do only if this edge flux is needed to update owned cells
         !
         if (cell1 <= nCells .or. cell2 <= nCells) then
 
            cell_list(1) = cell1
            cell_list(2) = cell2
            n = 2 
  
          !  add cells surrounding cell 1.  n is number of cells currently in list
            do i = 1, nEdgesOnCell(cell1)
               if (cellsOnCell(i,cell1) /= cell2) then
                  n = n + 1
                  cell_list(n) = cellsOnCell(i,cell1)
               end if
            end do
  
          !  add cells surrounding cell 2 (brute force approach)
            do iCell = 1, nEdgesOnCell(cell2)
               addcell = .true.
               do i=1,n
                  if (cell_list(i) == cellsOnCell(iCell,cell2)) addcell = .false.
               end do
               if (addcell) then
                  n = n+1
                  cell_list(n) = cellsOnCell(iCell,cell2)
               end if
            end do
  
  
            nAdvCellsForEdge(iEdge) = n
            do iCell = 1, nAdvCellsForEdge(iEdge)
               advCellsForEdge(iCell,iEdge) = cell_list(iCell)
            end do
  
          ! we have the ordered list, now construct coefficients
  
            adv_coefs(:,iEdge) = 0.
            adv_coefs_3rd(:,iEdge) = 0.
          
          ! pull together third and fourth order contributions to the flux
          ! first from cell1
  
            j_in = 0
            do j=1, n
               if( cell_list(j) == cell1 ) j_in = j
            end do
            adv_coefs    (j_in,iEdge) = adv_coefs    (j_in,iEdge) + deriv_two(1,1,iEdge)
            adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(1,1,iEdge)
  
            do iCell = 1, nEdgesOnCell(cell1)
               j_in = 0
               do j=1, n
                 if( cell_list(j) == cellsOnCell(iCell,cell1) ) j_in = j
               end do
               adv_coefs    (j_in,iEdge) = adv_coefs    (j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
               adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
            end do
  
          ! pull together third and fourth order contributions to the flux
          ! now from cell2
  
            j_in = 0
            do j=1, n
               if( cell_list(j) == cell2 ) j_in = j
            end do
            adv_coefs    (j_in,iEdge) = adv_coefs    (j_in,iEdge) + deriv_two(1,2,iEdge)
            adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(1,2,iEdge)
  
            do iCell = 1, nEdgesOnCell(cell2)
               j_in = 0
               do j=1, n
                  if( cell_list(j) == cellsOnCell(iCell,cell2) ) j_in = j
               end do
               adv_coefs    (j_in,iEdge) = adv_coefs    (j_in,iEdge) + deriv_two(iCell+1,2,iEdge)
               adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(iCell+1,2,iEdge)
            end do
  
            do j = 1,n
               adv_coefs    (j,iEdge) = - (dcEdge(iEdge) **2) * adv_coefs    (j,iEdge) / 12.
               adv_coefs_3rd(j,iEdge) = - (dcEdge(iEdge) **2) * adv_coefs_3rd(j,iEdge) / 12.
            end do
  
          ! 2nd order centered contribution - place this in the main flux weights
  
            j_in = 0
            do j=1, n
               if( cell_list(j) == cell1 ) j_in = j
            end do
            adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
  
            j_in = 0
            do j=1, n
               if( cell_list(j) == cell2 ) j_in = j
            end do
            adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
  
          !  multiply by edge length - thus the flux is just dt*ru times the results of the vector-vector multiply
  
            do j=1,n
               adv_coefs    (j,iEdge) = dvEdge(iEdge) * adv_coefs    (j,iEdge)
               adv_coefs_3rd(j,iEdge) = dvEdge(iEdge) * adv_coefs_3rd(j,iEdge)
            end do
 
         end if  ! only do for edges of owned-cells
         
      end do ! end loop over edges

   end subroutine atm_adv_coef_compression


   subroutine atm_couple_coef_3rd_order(mesh, configs)

      implicit none

      type (mpas_pool_type), intent(inout) :: mesh
      type (mpas_pool_type), intent(in) :: configs

      real (kind=RKIND), dimension(:,:), pointer :: adv_coefs_3rd
      real (kind=RKIND), dimension(:,:,:), pointer :: zb3_cell
      real (kind=RKIND), pointer :: config_coef_3rd_order

      call mpas_pool_get_array(mesh, 'zb3_cell', zb3_cell)
      call mpas_pool_get_array(mesh, 'adv_coefs_3rd', adv_coefs_3rd)

      call mpas_pool_get_config(configs, 'config_coef_3rd_order', config_coef_3rd_order)

      adv_coefs_3rd(:,:) = config_coef_3rd_order * adv_coefs_3rd(:,:)
      zb3_cell(:,:,:) = config_coef_3rd_order * zb3_cell(:,:,:)

   end subroutine atm_couple_coef_3rd_order


   !-----------------------------------------------------------------------
   !  routine mpas_atm_run_compatibility
   !
   !> \brief Checks input fields and options for compatibility
   !> \author Michael Duda
   !> \date   22 October 2018
   !> \details
   !>  This routine checks the input fields and run-time options provided
   !>  by the user for compatibility. For example, two run-time options may
   !>  be mutually exclusive, or an option may require that a certain input
   !>  field is provided.
   !>
   !>  A value of 0 is returned if there are no incompatibilities among
   !>  the provided input fields and run-time options, and a non-zero value
   !>  otherwise.
   !
   !-----------------------------------------------------------------------
   subroutine mpas_atm_run_compatibility(dminfo, blockList, streamManager, ierr)

#ifdef DO_PHYSICS
      use mpas_atmphys_control, only : physics_compatibility_check
#endif
      use mpas_atm_boundaries, only : mpas_atm_bdy_checks
      use atm_time_integration, only : mpas_atm_dynamics_checks

      implicit none

      type (dm_info), pointer :: dminfo
      type (block_type), pointer :: blockList
      type (MPAS_streamManager_type), pointer :: streamManager
      integer, intent(out) :: ierr

      integer :: local_ierr

      ierr = 0

#ifdef DO_PHYSICS
      !
      ! Physics specific checks found in physics/mpas_atmphys_control.F
      !
      call physics_compatibility_check(dminfo, blockList, streamManager, local_ierr)
      ierr = ierr + local_ierr
#endif

      !
      ! Checks for limited-area simulations
      !
      call mpas_atm_bdy_checks(dminfo, blockList, streamManager, local_ierr)
      ierr = ierr + local_ierr

      !
      ! Checks for dynamics options
      !
      call mpas_atm_dynamics_checks(dminfo, blockList, streamManager, local_ierr)
      ierr = ierr + local_ierr

   end subroutine mpas_atm_run_compatibility

end module atm_core

